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ABSTRACT 

The  report  considers  the  problem  of  target  motion  analysis  from  range  and 
range-rate  target  measurements.  The  motivation  for  this  work  comes  from  the 
need  to  track  the  target  motion  with  the  Ingara  Multi-Mode  Radar  during  an 
extended  data  collection  in  the  ISAR  mode.  The  report  makes  three  main 
contributions.  First,  the  theoretical  Cramer-Rao  bound  for  the  performance  of 
an  unbiased  range-only  tracking  algorithm  is  derived.  Second,  three  algorithms 
for  target  motion  analysis  (using  range  and  range-rate  measurements  only)  are 
developed  and  compared  to  the  theoretical  bounds  of  performance.  The  three 
algorithms  are:  the  Maximum  Likelihood  estimator,  the  Extended  Kalman 
filter  and  the  Regularised  Particle  filter.  Finally,  the  report  presents  the  results 
of  the  application  of  the  developed  theory  to  the  ISAR  data  collected  in  the 
recent  trials  with  the  Ingara  radar. 
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Target  Motion  Analysis  Using  Range-Only  Measurements: 
Algorithms,  Performance  and  Application  to  Ingara  ISAR 

Data 


EXECUTIVE  SUMMARY 

The  report  investigates  the  problem  of  target  motion  analysis  from  range  and  range- 
rate  measurements.  The  problem  is  of  importance  for  tracking  surface  vessels  while  col¬ 
lecting  data  in  the  ISAR  mode  of  the  Ingara  multi-mode  radar.  The  report  derives  the 
bounds  of  performance  and  considers  three  algorithms  for  target  motion  analysis.  The 
bounds  indicate  the  criterion  for  target  state  observability,  which  can  be  summarised  as  a 
requirement  for  the  observer  to  “outmanoeuvre”  the  target.  The  proposed  algorithms  for 
target  motion  analysis  include  the  Maximum  Likelihood  estimator,  the  Extended  Kalman 
filter  (EKF)  with  and  without  angle  parametrisation  and  the  Regularised  Particle  filter. 
Considering  both  the  statistical  and  computational  performance  of  the  proposed  algo¬ 
rithms  and  based  on  extensive  performance  evaluation  with  simulated  data,  it  appears 
that  the  angle-parametrised  EKF  is  the  preferred  choice  for  implementation  in  an  opera¬ 
tional  system. 

The  discussed  algorithms  for  target  motion  analysis  have  been  applied  successfully 
to  a  set  of  real  ISAR  data  (range  profiles)  collected  in  the  recent  trials  with  the  Ingara 
radar.  The  results  indicate  that  target  motion  analysis  based  on  range  and  range-rate 
measurements  can  converge  towards  a  steady  state  in  less  than  three  minutes  for  a  typical 
scenario  of  interest. 
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1  Introduction 


The  problem  of  target  tracking  or  target  motion  analysis  (TMA)  using  range  and  range- 
rate  measurements  is  motivated  by  its  application  to  maritime  surveillance  with  the  Ingara 
Multi-Mode  Radar  [5].  In  particular,  one  of  the  maritime  surveillance  modes  of  Ingara  is 
the  Inverse  Synthetic  Aperture  Radar  (ISAR)  mode  [7],  in  which  the  radar  collects  high- 
resolution  range  profiles  of  the  target,  processed  in  real-time  to  generate  ISAR  images.  The 
problem,  however,  is  that  during  the  collection  of  the  range-profiles,  the  antenna  spotlights 
the  location  on  the  sea-surface  where  the  target  was  initially  detected,  while  in  reality  the 
target  is  moving  and  at  some  stage  disappears  from  the  radar  antenna  beam  in  azimuth 
or  the  sampled  swath  in  range.  Therefore  there  is  a  need  to  track  the  target  in  the  ISAR 
mode  using  only  the  information  contained  in  the  ISAR  range  profiles.  This  information 
can  be  presented  in  the  form  of  the  target  range  and  range-rate  measurements  for  tracking. 
If  target  tracking  is  possible  with  these  measurements,  the  output  of  the  tracker  could  be 
used  to  automatically  control  the  range  gate  setting  and  antenna  pointing. 

Tracking  targets  with  range  only  measurements  has  attracted  very  little  research  in¬ 
terest  in  the  past.  To  our  best  knowledge,  the  only  publication  devoted  to  this  problem, 
[19],  discusses  the  conditions  for  target  observability  from  range-only  measurements.  It 
concludes  that  if  the  target  is  moving  at  a  constant  acceleration  the  observer  must  be 
moving  with  a  non-zero  jerk  in  order  to  observe  the  target.  This  condition  appears  to  be 
identical  to  the  observability  criterion  for  the  related  and  extensively  studied  problem  of 
angle-only  target  motion  analysis  [10],  [4,  Ch.5]. 

This  report  makes  three  main  contributions.  First  it  derives  the  theoretical  Cramer- 
Rao  (CR)  bounds  of  performance  for  the  target  tracking  problem  with  range  and  range-rate 
measurements.  The  bounds,  which  define  the  best  achievable  performance,  depend  on  the 
target-observer  geometry,  the  accuracy  of  measurements  and  the  sampling  interval.  The 
derivation  assumes  unity  probability  of  detection  and  zero  probability  of  false  alarm  (i.e. 
no  uncertainty  with  respect  to  the  origin  of  measurements)  as  well  as  non-manoeuvering 
target  motion.  The  bounds  are  derived  for  the  recursive  type  estimators  with  prior  infor¬ 
mation  and  automatically  provide  the  criterion  for  observability.  The  second  contribution 
of  the  report  is  the  design  of  three  algorithms  for  state  estimation  and  their  compari¬ 
son  with  the  derived  CR  bounds.  The  problem  of  range-only  tracking  is  a  non-linear 
problem,  and  the  proposed  algorithms  include:  (i)  the  Maximum  likelihood  estimator 
(MLE)  over  a  cumulative  measurement  set;  (ii)  the  Extended  Kalman  filter  (EKF)  with 
angle-parametrisation  and  (iii)  the  Regularised  Particle  filter  (RPF).  The  MLE  is  a  batch 
algorithm,  but  in  this  application  it  is  applied  to  a  set  of  measurements  collected  within 
a  sliding  time  interval.  Both  the  EKF  and  the  RPF  are  recursive  algorithms.  The  third 
contribution  of  the  report  is  the  application  of  the  developed  range-only  tracking  theory 
to  a  set  of  real  ISAR  data  collected  in  the  recent  trials  with  the  Ingara  radar. 

The  report  is  organised  as  follows.  Section  2  describes  the  motivation  for  range-only 
tracking  in  ISAR  imaging  and  provides  a  mathematical  formulation  of  the  problem.  The 
CR  lower  bounds  of  performance  are  derived  and  discussed  in  Section  3.  Section  4  presents 
the  three  proposed  tracking  algorithms.  Their  statistical  and  computational  performance 
is  compared  in  Section  5.  Section  6  is  devoted  to  the  application  of  the  range-only  TMA 
algorithms  to  a  set  of  ISAR  range-profiles  collected  by  the  Ingara  radar.  Section  7  sum- 
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marises  the  main  results  of  the  report. 


2  Problem  description 

2.1  Background 

Two  of  the  maritime  surveillance  modes  of  the  Ingara  radar  are  the  Scan  mode  [15] 
and  the  ISAR  mode  [7].  In  typical  wide  area  maritime  surveillance,  an  airborne  observer 
is  using  the  radar  Scan  mode  to  detect  and  locate  maritime  surface  vessels.  Once  the 
location  of  a  target  is  known,  the  radar  operator  switches  the  radar  into  the  ISAR  mode 
to  collect  high  resolution  range  profiles  of  the  target.  The  range  profiles  are  processed 
into  the  ISAR  images  and  shown  to  the  operator  in  real-time.  In  the  ISAR  mode,  the 
radar  antenna  is  controlled  in  such  a  manner  that  it  “spotlights”  the  ocean  surface  at 
the  target  coordinates  which  were  obtained  from  the  initial  detection  in  the  radar  Scan 
mode.  If  the  target  is  moving,  however,  at  some  point  of  time  it  will  move  out  of  the 
radar  antenna  beam  either  in  the  cross-range  direction  or  the  sampled  swath  in  the  range 
direction.  Because  we  wish  to  continuously  collect  the  data  in  the  ISAR  mode,  we  would 
like  to  be  able  to  track  the  target  using  only  the  information  contained  in  the  ISAR  range 
profiles.  This  would  enable  the  radar  antenna  beam  to  be  steered  always  towards  the 
target.  The  information  contained  in  the  high-resolution  range  profiles  can  be  converted 
to  the  target  range  and  range-rate  measurements  by  linear  regression  of  the  target  range 
centroids,  applied  over  a  sliding  time  window. 

This  motivates  the  theoretical  consideration  of  the  problem  of  target  motion  analysis 
or  target  tracking,  using  only  the  available  measurements  of  target  range  and  range-rate. 
Before  developing  the  tracking  algorithms  the  fundamental  questions  to  be  answered  are: 
is  the  range  tracking  feasible  (the  observability  criterion)  and  how  fast  the  algorithms  can 
converge.  Note  that  for  the  ISAR  data  collection  purposes,  the  observer  is  usually  flying 
along  a  circular  trajectory,  as  illustrated  in  Figure  1.  This  is  the  type  of  the  observer 
trajectory  we  will  consider  in  the  report. 
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2.2  Mathematical  formulation 

The  target  is  assumed  to  move  at  a  constant  velocity  along  a  straight  line  in  the  x-y 
plane  (Figure  1).  Although  in  reality  the  target-observer  scenario  is  three-dimensional,  it 
suffices  to  define  and  study  the  problem  in  two  dimensions,  because  the  slant  range  and 
range-rate  can  be  readily  converted  to  the  ground  range  and  range-rate1.  For  the  typical 
ranges  of  interest  the  curvature  of  the  Earth  can  be  neglected  and  hence  we  assume  the 
local  tangential  plane  (flat  Earth)  model. 

The  basic  problem  is  to  estimate  the  target  kinematic  state  (position,  heading,  speed) 
from  noise  corrupted  measurements.  The  target  kinematic  state  can  be  fully  described  by 
the  state  vector  defined  in  the  discrete-time  as: 

4  =  [4  4  4  vif  (i) 

where  T  denotes  matrix  transpose,  x\  and  y\  are  the  Cartesian  target  coordinates  at  time 
index  k  and  xj.  and  yk  are  their  respective  derivatives  (velocities).  For  the  assumed  target 
motion  model,  xk  and  yk  are  constant  in  time.  Knowing  the  target  state  vector  it  is  then 
straightforward  to  calculate  the  target  heading  or  speed.  The  observer  state,  defined  as 

4  =  [4  4  vl  y°k}T  (2) 

is  assumed  to  be  known  (information  supplied  by  an  on-board  Inertial  Navigation  System). 
The  relative  state  vector  is  given  by: 

sk  =  stk-s°k  =  [xk  Xk  yk  yk]T.  (3) 

The  available  measurements  of  target  range  and  range-rate  are  available  at  time  in¬ 
stants 


tk  =  t0  +  Td  4-  (k  -  1  )TS,  A;  =  1,2, ...  (4) 

where  Ts  is  the  sampling  interval  in  the  ISAR  mode  and  Td  »  Ts  is  the  period  of  time 
required  for  the  radar  to  switch  from  the  Scan  mode  to  the  ISAR  mode.  The  assumption 
is  that  the  measurement  origin  is  unambiguous2  and  that  the  target  is  detected  with  the 
probability  1.  The  measurement  vector  at  time  instants  tk  is  defined  by: 


z k  = 


rk 

rk 


(*  =  1,2,...) 


and  the  measurement  equation  is  given  by: 


(5) 


zk  =  h'{sk)  +  wk 


"  rk  “ 

hr  (sk) 

.  ^k  . 

hy{sk) 

(6) 

(7) 


The  transformation  from  slant  to  ground  range  and  range-rate  is  straighforward  since  the  target  is 
always  on  the  surface  of  the  Earth  and  the  observer  height  is  constant  and  known  during  the  observation 
interval. 

2 Measurement  origin  ambiguity  could  be  caused  for  example  by  multiple  targets  or  possibly  false  mea¬ 
surements. 
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where 


hr  (sfc) 
hr  (s/c) 


=  y/4  +  Vk 

Xk  Xk  +  yk  Vk 


\[x 


+  vl 


(8) 

(9) 


The  noise  process  w*.  in  (6)  is  i.i.d.  zero-mean  Gaussian  with  covariance  R,  that  is 
wk  ~  «A/*(0;R). 

The  information  about  the  target  state  at  time  to  is  available  in  form  of  the  measure¬ 
ment  of  the  target  location  (range  and  azimuth)  provided  by  the  radar  in  the  Scan  mode. 
This  measurement  is  denoted  by 


zo 


ro 

do 


and  the  measurement  equation  at  to  is  given  by: 

z0  =  h"(s0)  +  w0 
"o  _  hr{so)  , 

>oJ  -  Im«)J+  0 

where  hT{ )  was  defined  by  (8)  and 


he{ so)  =  arctan  (  —  )  . 

V  3/o  y 


(10) 

(11) 

(12) 

(13) 


In  (13),  xq  and  yo  are  the  relative  target  Cartesian  coordinates  at  k  =  0  and  can  be  worked 
out  from  the  measurement  at  to-,  zo  (which  is  in  polar  coordinates).  Note  that  azimuth  0o 
is  defined  as  the  angle  from  the  observer  to  the  target,  referenced  clockwise  positive  to  the 
y-axis.  The  noise  process  wo  in  (11)  is  independent  from  w k  =  1,2,...  and  distributed 
as  w0  ~  J\f( 0;  Ro). 

The  state  equation  for  this  problem  can  be  written  as: 

1  “  ~  Ufc+i,* 

where  is  the  transition  matrix  defined  as 

1  Tk  0  0 

0  10  0 
0  0  1  Tk 

0  0  0  1 


(14) 


(15) 


Tk  is  the  sampling  interval,  given  by 

t‘={ 

and 


Td 

Ts 


k  =  0 

k  =  1,2,... 


(16) 


U 


k+l,k  — 


*2+1  -  *2  -  Tk'X 


*2+i  -  xl 


TkVk 


Vk+ 1  -  Vk 

y°k+ 1  -  v°k 

is  a  vector  of  deterministic  inputs  which  account  for  the  observer  acceleration.  Note  that 
the  model  of  the  state  equation  (14)  is  purely  deterministic. 
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3  Cramer-Rao  Lower  Bounds 

3.1  Derivation  of  the  Bounds 

The  Cramer-Rao  lower  bound  (CRLB)  for  an  unbiased  estimator  s£  of  the  target  state 
vector  at  time  tk,  is  given  by  the  inverse  of  the  Fisher  information  matrix  (FIM), 

e{(4-4)  (4  -  4)r}  >  J*1.  (*  =  1,2,...).  (17) 

Since  the  observer  motion  is  deterministic,  the  same  CRLB  as  in  (17)  applies  to  an  unbiased 
estimator  of  the  relative  target  state  vector  sk.  The  FIM  in  (17)  is  defined  as  [20]: 

J*  =  E{[vA(s*)][vA(s£)f}  (18) 

where  A(sj[)  is  the  negative  log-likelihood  function 

A(st)  =  -logp(Zfc|s£),  (19) 

and  V  is  the  gradient  with  respect  to  the  target  state  vector  s[.  The  term  Zk  in  (19) 
represents  the  collection  of  all  measurements  up  to  time  tfC)  i.e.  Z&  —  {z i}f=o- 

Due  to  the  assumption  we  made  about  the  independence  of  measurement  noise,  the 
likelihood  function  can  be  expressed  by  a  product: 

k 

p(zfcl4)  =  rb(z*ls*)-  (2°) 

i=0 

Furthermore,  since  the  measurement  noise  is  zero-mean  Gaussian 


'V(h"[s0(4)];Ro)  i  =  0 

P(*i|8fc)  =  * 

(21) 

,V(h'[Si(s^.)];R)  *=1,2,...,  A 

27TV/jR^eXP  [  z(Z0  “  11  [So(8*)])  R-0  (z0  -  h  [s0(sfc)])] 

i  =  0 

=  < 

(22) 

1  2^  exp  [  -  h'[Si(st)])r  R  (z,-  -  h'[Si(s‘ )])] 

i  =  1,2, . ..  ,fc. 

Functions  h'(  )  and  h"( )  in  (22)  were  defined  in  Section  2.2  as 

with  hr(  ),  hf{  )  and  h$(  )  given  by  eqs.(8),  (9)  and  (13)  respectively.  The  crucial  point 

here  is  to  observe  that  the  arguments  of  h'(  )  and  h"(  )  have  to  be  expressed  as  the 

respective  functions  of  s£.  This  means  that 

1.  xo  and  yo  in  (13)  should  be  expressed  as 

*o(s£)  =  x\  -  {tk  -  t0)x{  -  x°0  (24) 

W(4)  =  Vk  ~  (h  ~  to)y{  -  y°0  (25) 

which  follows  from  the  earlier  assumption  that  ±1  —  const  and  y\  -  const.  Note 
from  (4)  that  tk  -  t0  =  Td  -f  (k  -  1  )TS. 
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2.  Xi,  yi,  and  yi  (i  =  1,2,...,  k)  which  according  to  (22)  are  to  be  used  in  eqs.(8) 
and  (9)  have  to  be  expressed  as: 


*<(4) 

=  x{-  (k-i)  Ts  x{  -  x\ 

(26) 

Vi(  4) 

=  y{-(k-  *)  Ts  y{  -  yf 

(27) 

M  4) 

=  **“*? 

(28) 

w(4) 

=  y{-y°i- 

(29) 

Again  eqs.(26)-(29)  follow  from  x\  =  const  and  yk  =  const. 

From  (19)  and  (22)  it  follows  that  the  log-likelihood  function  can  be  expressed  as: 

A(4)  =  C+  \  {®o  -  h"[s0(s|)]}  Rq1  {z0  -  h"[s0(s^)]}T  + 

1  ^ 

2  E  {*  -  h'M4)]}  {*<  -  h'[Si(4)]}T  (30) 

Z— 1 

k 

=  c + Ao  +  y  ^  A  i 

2=1 

where  C  is  a  constant.  Next  we  need  to  find  the  gradient  of  the  log-likelihood  function  A 
expressed  by  (30),  which  using  the  basic  rules  of  matrix  algebra  [12]  yields: 


k 

VA(s£)  =  VA0  +  ^VA*  (31) 

2=1 

with 

VA0  =  -V[h"]T  Rq  1  [z0  -  h"]  (32) 

VAi  =  -V[h']T  R-1  [zj  -  h"]  (*  =  1,2,...,*)  (33) 

The  FIM  of  eq.(17)  is  then 

{•  k  I  r  it 

^VAi  X^VAi  l  (34) 

.2=0  _  .2=0  _  J 

which  due  to  the  mutual  independence  of  measurements  leads  to 

J„  =  [V[h"f ]  R,'  [V[h”f]3'  +  £  [v[hf]  R-1  [v[hff  (35) 

2  =  1 


It  remains  to  find  the  analytic  expressions  for 
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Using  the  differentiation  rules,  the  unknown  terms  in  (36)  can  be  shown  to  be: 


dhr[si{sl)} 

Xi 

dxk 

s/xi+Vi 

dhr[Si(  sjj.)] 

dii 

= 

-t,k-  i)r,aft^s‘f4)1 

dhrlsiisj.)} 

Vi 

H 

'jA  +  y’i 

dhr  [sj(s^.)] 

= 

-lk-i)T.dh' ff**1 

dh^  [si  (sjj,)] 

xi\Jxf  +  Vi  (xiXi  +  yiVi)  7=^ 
v  vxi+Vi 

dxk 

xf  +  Vi 

dhf[si{sl)] 

[xi  (k  i)TsXi]  ^ xf  +  yf  +  (Xi  Xi  +  yiyi )  ^  ^ | 

dil 

xf  +  yf 

5hr[st(s^.)] 

yi\jx2i  +  Vi  +  Villi )  v/J’+y2 

dVk 

x?  +  yf 

3MM4)] 

[Vi  ( k  l)Tsih]  J xf  +  yf  +  ( Xi  Xi  +  yi  y,) 

dyi 

xf  +  yf 

dhr[ s0(s^)] 

x0 

dx{ 

Vxo  +  Vo 

dhT  [so  (sjj.)] 

dik 

= 

dhr[sQ(sfk)] 

yo 

dyl 

Vxf+vl 

dhr[so{sl)] 

dVk 

= 

dyi 

dhe[  so(s|.)] 

yo 

dxk 

Xq  +  Vo 

9he  [s0  (a*k)] 
d±\ 

= 

-[(*-  Dr.  +  r/^W)] 

dheisoisi)] 

^0 

dVk 

xf  +  yf 

dhe{  so(s^)) 
dyi 

= 

dyi 

(37) 

(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 


The  FIM  given  by  (34)  depends  on:  (i)  the  geometry  of  the  considered  scenario  (target  and 
observer  trajectories);  (ii)  measurements  accuracy  (defined  by  the  covariance  matrices  Ro 
and  R);  (iii)  time  delay  T'd  and  the  sampling  interval  Ts.  The  CRLBs  of  the  state  vector 
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4  are  calculated  as  the  diagonal  elements  of  the  inverse  of  the  FIM,  i.e. 

CRLB{  [4].}  =  [J,-1] . .  (j  =  1, 2, 3,4),  (53) 

where  [sj(]  .  refers  to  the  j-th  component  of  vector  s(..  Note  that  the  target  observability 
implies  that  the  FIM  has  an  inverse  matrix. 


3.2  Analysis  of  the  Bounds 

Consider  a  scenario  shown  in  Figure  2.  The  target  is  moving  towards  the  north-east 
with  a  speed  of  8  m/s.  The  observer  is  flying  at  a  speed  of  150  m/s  along  a  circular 
path  with  a  radius  of  approximately  15  km.  At  time  to ,  the  observer  is  at  a  location 
with  coordinates  (50km,  35km)  in  Figure  2.  The  covariances  are  assumed  to  be  R0  — 
diag  [400m2,  (30)2]  and  R  =  diag  [400m2,  4m2/s2],  with  a  time  delay  of  Td  =  40s  and 
sampling  interval  Ts  =  Is.  The  CRLB’s  (square-rooted)  of  unbiased  estimators  xk,  Vk > 


Figure  2:  A  tracking  scenario  used  for  the  analysis  of  Cramer-Rao  lower  bounds 

xk  and  yk  (k  =  1, 2, . . . ),  are  shown  in  Figure  3  (starting  from  time  instant  t\  =  t0  +  Td). 
Note  that  the  initial  range  estimate  is  fairly  accurate  while  the  initial  angular  estimate 
is  fairly  poor.  This  translates  in  this  particular  observer-target  scenario,  into  a  more 
accurate  initial  estimate  of  the  target  position  in  y  than  along  the  x  axis.  With  estimates 
of  velocities  x  and  y  the  situation  is  opposite. 

Next  we  analyse  the  effect  of  the  following  parameters  on  the  CRLBs:  the  time  delay 
between  the  radar  modes  (Td),  the  sampling  interval  (T5),  the  radius  of  the  observer 
circular  path  p  and  the  standard  deviation  of  the  range-rate  measurements  =  y^R^  • 
The  results  are  shown  in  Figure  4  for  the  same  scenario  and  the  parameters  used  in  Figures 
2  and  3,  except  for  those  parameters  whose  effect  on  the  CRLBs  is  being  analysed. 

Somewhat  surprisingly,  the  CRLBs  in  Figure  4.  (a)  and  (b)  suggest  that  longer  radar 
mode  switch  delays,  Td,  speed  up  the  convergence  of  the  range-only  tracking.  The  expla¬ 
nation  is  that  by  extending  Td  we  effectively  increase  the  duration  of  the  total  observation 
period  tk  —  to,  which  facilitates  the  state  estimation. 
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Figure  3:  Cramer-Rao  lower  bounds:  (a)  position  Xk  (solid  line)  and  y *  (dashed  line);  (b) 
velocity  x (solid  line)  and  (dashed  line) 


By  reducing  the  sampling  interval,  the  total  number  of  measurements  is  increased,  and 
the  estimates  of  the  target  position  in  the  x  and  y  directions  are  expected  to  be  better, 
as  confirmed  by  the  CRLBs  in  Figure  4.(c)  and  (d).  Note  that  this  would  increase  the 
computational  load  of  the  algorithms. 

The  effect  of  varying  the  radius  of  the  circular  observer  path,  p,  is  shown  in  Figure 
4.(e)  and  (f).  Since  the  target  is  moving  with  a  constant  velocity  along  a  straight  line, 
the  observability  criterion  [19]  requires  the  observer  to  have  a  non-zero  acceleration.  The 
smaller  the  radius  of  the  observer  circular  path,  the  larger  the  component  of  the  observer 
acceleration,  and  consequently  the  faster  the  convergence  of  the  state  estimates  from  range 
and  range-rate  only  measurements.  It  has  been  verified  that  if  the  observer  is  moving  along 
a  straight  line  with  a  constant  velocity  (same  as  the  target),  the  CRLBs  are  monotonically 
increasing,  thus  confirming  the  observability  requirement  for  non-zero  acceleration. 

The  significance  of  the  range-rate  measurements,  r&,  is  investigated  in  Figure  4.(g)  and 
(h).  Two  cases  are  considered,  the  first  with  totally  imprecise  range-rate  measurements 
(i.e.  very  large  [R]2,2,  mimicking  the  absence  of  r *.)  and  the  second  with  precise  range-rate 
measurements  Gf  —  1  m/s.  The  CR  bounds  indicate  that  by  using  the  precise  range-rate 
measurements,  the  algorithm  convergence  can  be  improved.  Since  the  algorithm  conver¬ 
gence  is  of  crucial  importance  in  the  application  of  interest,  the  range-rate  measurements 
will  be  used  for  target  tracking. 


3.3  Bounds  with  the  use  of  prior  knowledge 

The  recursive  Bayesian  estimators,  such  as  the  Extended  Kalman  filter  or  the  particle 
filter  (described  in  depth  in  Section  4),  in  addition  to  all  available  measurements,  Z*.,  make 
use  of  prior  knowledge  about  the  target  motion  for  tracking.  For  the  problem  considered 
in  this  report,  this  could  be  prior  knowledge  of  the  maximum  target  speed. 

The  recursive  Bayesian  estimators  typically  estimate  the  relative  target  state  vector 
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and  initialise  the  target  track  using  the  first  measurement  and  the  prior  knowledge.  For 
the  problem  considered  in  this  report  the  initial  measurement  at  to,  zo  =  [ro,  0o]T,  is  in 
polar  coordinates.  The  initial  target  state  and  its  covariance  are  calculated  by  conversion 
from  the  polar  to  the  Cartesian  coordinates3,  assuming  that  the  covariance  matrix  R0  is 
diagonal  with  elements  aT  =  -\/[Ro]i,i  and  —  ^/[Ro]2,2-  The  initial  state  vector  is  then 

So|o  =  sin  0O  -Ag  ro  cos  60  -  (54) 

where  Xq  and  are  the  initial  observer  velocity  components  defined  in  eq.(2).  The  co- 
variance  matrix  of  the  initial  estimate  is  given  by 

Pn  0  P\3  0 

p  ,  _  0  P22  0  0 

0,0  Psi  o  P33  0  •  <55) 

0  0  0  P44  _ 

Note  the  initial  target  velocity  components  in  (54)  are  set  to  zero,  while  their  variance  is 
given  by  the  elements  P22  and  P44  of  P0|0  in  (55),  i.e. 

-P22  =  P44  =  cr2.  (56) 

The  remaining  elements  of  P0|0  are  as  follows  [8,  p.155]: 

P11  =  ro<70  cos2  0q  +  <r2  sin2  60 
P33  =  rgof  sin2  <90  +  o'2  cos2  0q 
P13  —  P31  =  {&r  ~  ro(7o)  sin0o  cos  00- 

The  prior  target  density  at  to  is  assumed  to  be  Gaussian,  that  is  so  ~  A/"(so|o, Po|o)-  The 
contribution  of  this  prior  density  now  has  to  be  worked  out  at  the  time  index  k.  Due  to 
the  purely  deterministic  target  motion,  the  following  relationship  holds: 


where 


sk  =  Ffc  s0 


1  (k-l)Ts  +  Td  0  0 

0  10  0 

0  0  1  {k-l)T,+Td 

0  0  0  1 


Then  it  follows  (see  for  instance  [13,  Theorem  2.11])  that 

s*~^(FfcSo|0>  FfeP0|o[F*f). 

The  FIM  at  time  index  k  for  the  recursive  Bayesian  estimators  can  be  expressed  as 
[20,  pg.84]  : 


3The  so  called  “debiased  consistent”  conversion  [14]  is  recommended  for  large  values  of  <j$. 
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where  [20,  pg.85] 


Jfc  =  [FftP0|0[Ffc]r' 


(59) 


is  due  to  the  initial  measurement  and  the  prior  knowledge.  The  component  of  the  FIM 
due  to  measurements,  Jj^,  follows  from  (35): 


Jf  =  £[v[h']T]  R-1  [v[h']Tf. 

i= 1 


(60) 


Figure  5  compares  the  CRLB  with  and  without  the  use  of  prior  knowledge  ( x  position 
only).  The  trajectories,  measurement  covariances,  sampling  interval  and  the  initial  delay 
are  the  same  as  for  the  CRLBs  shown  in  Figure  3.  The  CRLBs  in  Figure  5  are  shown 
for  values  of  the  standard  deviation  of  velocity  (av  in  eq.(56))  of  5  m/s,  10  m/s  and  20 
m/s.  Observe  from  Figure  5  that  the  CRLB  which  uses  prior  knowledge  is  initially  lower, 


Figure  5:  Comparison  of  square-rooted  Cramer-Rao  lower  bounds  ( x  position  only)  with  and 
without  the  use  of  prior  information 

although  the  significance  of  prior  information  diminishes  over  time  and  the  bounds  become 
identical  for  k  >  100.  The  quality  of  prior  information  is  measured  by  its  variance.  Note 
that  by  increasing  the  value  of  av  (i.e.  by  reducing  the  quality  of  prior  information),  the 
CRLB  with  prior  knowledge  approaches  the  one  with  no  prior  knowledge.  In  summary, 
for  large  values  of  av ,  the  CRLB  calculated  using  the  FIM  of  eq.(35)  and  the  CRLB 
calculated  using  the  FIM  of  eq.(58)  become  identical.  Note  that  the  practical  value  of  av 
for  range-only  tracking  of  surface  vessels  is  av  ^  10  m/s  (dotted  line  in  Figure  5)  where 
prior  knowledge  makes  a  measurable  impact  on  the  initial  performance. 
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4  Tracking  Algorithms 

4.1  Maximum  Likelihood  Estimation 

The  Maximum  Likelihood  Estimator  (MLE)  [20]  of  the  target  state  vector  s£  ==[a;jj., 
ijj.,  yk]T  at  time  tk  is  defined  as  the  vector  s£  that  maximises  the  likelihood  func¬ 
tion  p(Zfc|sjj.),  given  by  eqs.(20).  This  is  equivalent  to  minimisation  of  the  negative  log- 
likelihood  function  A(s£),  given  by  (30).  The  implementation  is  easier  with  the  latter 
approach,  hence  we  find  the  MLE  as 

4, mle  =  arg  min  A(s£).  (61) 

Sk 

Minimisation  is,  in  general,  performed  by  numerical  methods  [18].  Our  implementation  is 
based  on  the  MATLAB©  built-in  routine  fmins.  Note  that  the  MLE  is  a  batch  algorithm 
-  it  operates  on  the  accumulated  set  of  all  previous  measurements.  For  on-line  data 
processing  applications  (such  as  the  one  considered  in  this  report),  the  MLE  can  be  applied 
over  a  sliding  window  of  accumulated  measurements.  From  the  theory  it  is  known  [20] 
that  the  MLE  is  an  efficient  estimator,  and  hence  its  variance  should  meet  the  theoretical 
CR  lower  bound  derived  in  Section  3.1  (no  use  of  prior  knowledge). 


4.2  Angle-Parametrised  Extended  Kalman  filter 

4.2.1  Extended  Kalman  filter 


Target  motion  analysis  using  range  and  range-rate  measurements  only,  is  a  non-linear 
dynamic  state  estimation  problem  because  the  measurement  equation  (6)  is  non-linear. 
The  linear  filtering  algorithms  such  as  the  Kalman  filter  are  therefore  inappropriate,  and 
the  conventional  approach  is  to  approximate  eq.(6)  by  a  series  expansion  and  then  to 
use  an  equivalent  measurement  matrix  in  the  ordinary  Kalman  filter  equations.  Such  an 
estimator,  known  as  the  Extended  Kalman  filter  (EKF)  [2],  is  developed  for  the  range- 
only  tracking  problem,  with  the  first  order  series  expansion  (linearisation)  of  the  non-linear 
measurement  equation. 

The  recursive  equations  of  the  EKF  are  presented  below.  They  describe  how  to  evaluate 
the  relative  state  estimate  s^+1|fc+1  at  tk+i  and  its  covariance  matrix  P^+i|A:+ii  given  the 
measurement  Zfc+i,  and  the  estimate  of  the  relative  state  vector  sk  at  tkl  with  its  covariance 
matrix  Pk\ k.  The  state  prediction  equation  follows  from  (14): 

®&4T|A;  —  F&  &k\k  ~  U&+1,A:  (62) 

while  the  covariance  prediction  is  given  by 

P/c+l|/c  =  ^kPk\k^k' 

The  measurement  prediction  is  given  by: 


£fc+i| k  = 


_  ^k-\-l\k  . 


=  h'(Sfe+l|fc) 
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with 


rk+l\k 


xk+i\k +  /4+ i|fe> 


(63) 


^k+l\k  ~ 


3'k+l\k%k+l\k  d~  V k-{-l\kV k-{-l\k 

fk+l\k 


The  Kalman  gain  matrix  can  be  evaluated  as: 


(64) 


Wfc+1  =  Pfc+i|fcHfc+i  [H*+iPA.+1|A.H^+1  +  R]  1  >  (65) 


where  R  is  the  covariance  matrix  of  measurement  noise  in  eq.(6),  and  H^+i  is  the  linearised 
measurement  matrix,  evaluated  at  the  predicted  state,  defined  by: 


d[h']i 

/in 

hi2 

^13 

hu 

d[s]j 

S—Sk+I\k 

/l21 

h  22 

^23 

/l24 

(66) 


The  elements  of  the  matrix  in  (66)  can  be  easily  computed  by  differentiation: 


J  &k+l\k 

mi  =  - - > 

rk+l\k 


hn  =  0,  h13  =  h14  =  0 

rfc+i|fc 


(67) 


^21 

^22 

^23 

h24 


_  ^'k-\-l\k^'k-\-l\k  3'k-\-l\k'i~k+l\k 


r 2 

f  k+l\k 


fe¬ 


ll 


_  rk+l\kyk+l\k  yk-\-l\k^k-\-l\k 


'  fc+l|fc 


=  fei3. 


The  updated  state  and  its  covariance  matrix  are  given  by 

Sfc+1|A;+1  =  &k+l\k  +  W^+i[z^+i  —  Zk+l\k\ 

PA;  +  l|/fc+l  =  [I  “  W&+lHfc+i]Pfc+i|fc 

where  I  is  the  4x4  identity  matrix. 


(68) 

(69) 

(70) 


4.2.2  Angle  Parametrisation 

If  the  initial  azimuth  measurement  is  characterised  by  a  fairly  large  variance,  the 
convergence  of  the  range-only  TMA  with  an  EKF  can  be  improved  by  running  in  parallel 
a  set  of  weighted  EKFs,  each  with  a  different  initial  azimuth  value.  This  type  of  a  filter 
in  the  report  will  be  referred  to  as  the  angle-parametrised  EKF  (AP-EKF).  The  idea  used 
in  the  development  of  the  AP-EKF  closely  resembles  that  of  the  range-parametrised  EKF 
for  angle-only  tracking  [17],  [1]. 

Recall  that  the  initial  azimuth  measurement  is  9q  with  standard  deviation  oq.  We  wish 
to  construct  Nf  Extended  Kalman  filters,  each  with  an  initial  angular  accuracy  Nf  times 
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better  than  that  of  a  single  EKF.  First  define  an  interval  in  the  angle  domain  with  the 
limits  6m in  =  9q  —  3<r q  and  0max  =  #o  +  3<j0.  The  interval  is  divided  into  Nf  subintervals  of 
width  A 6  =  {dmax-6min)/Nf.  The  initial  azimuth  angle  for  EKF  i  is  chosen  to  correspond 
the  mid-point  of  the  interval  z,  i.e. 

0i  =  ^min  +  (i  -  0.5)  A 0,  i  —  1, . . .  ,iV/, 


with  the  standard  deviation  set  to  cr#.  —  a#/Nf. 

The  next  step  is  to  calculate  the  weights  associated  with  each  EKF  at  time  £&+ 1-  This 
can  be  done  recursively  using  Bayes’  rule  as: 


i  _  p(gfc+il»)*4 

*+1  Ej2iP(zfc+ili)wfc 


(71) 


where  p(zjfc+i|z)  is  the  likelihood  of  measurement  zjt+i,  given  that  the  i- th  EKF  is  the 
correct  one.  Assuming  Gaussian  statistics,  the  likelihood  p(z*.+i|z)  can  be  computed  as: 


P(*fc+il*)  =  0'  i  {- \ (zfc+i  -  z*+i|Jfc)  (SUi)  1  (z*+i  “  zfc+i|fc)  }  (72) 

where  zj,+1|fc  is  the  predicted  measurement  vector  at  k  +  1  for  EKF  i  and  Sfc+1  is  the 
innovation  covariance  for  EKF  z,  given  by: 

Si+i  =  Hfc+1P*+1|fcH^+1  +  R.  (73) 

In  (73),  Hfc+i  is  the  linearised  measurement  vector  defined  by  eqs.(66)-(68),  pi+i\k is  the 
covariance  matrix  of  the  predicted  state  of  the  z-th  filter. 

The  initial  weights  wz0  are  computed  assuming  that  the  true  initial  azimuth  obeys  a 
Gaussian  distribution  with  mean  9q  and  variance  oj,  i.e. 


w0  = 


\p2/KGQ 


reR 

Je\ 


itzMl 
2ae  d6 


(74) 


where  9\  =  0min  +  (i  -  1)  A0  and  0^  =  0}.  +  A0  are  the  lower  and  upper  limits  of  the  z-th 
subinterval,  respectively. 

Suppose  the  updated  state  estimate  of  EKF  i  at  tk  is  denoted  by  s^.  The  com¬ 
bined  AP-EKF  state  vector  and  its  covariance  are  calculated  using  the  Gaussian  mixture 
formulae  [2,  p.47]: 


Nf 

Sk\k  =  YlWk§l\k 

i= 1 
Nf 

Wk  [Pfc|fc  +  fc|fc  ~  §*|&)(Sft|A.  -  Sk\k)']  ■ 


(75) 

(76) 


The  improved  tracking  performance  of  the  AP-EKF  tracker  is  achieved  by  using  Nf 
independent  trackers,  each  with  a  much  smaller  initial  angular  measurement  variance  than 
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that  of  a  single  EKF.  This  improvement  is  achieved  at  the  expense  of  an  ATy-fold  increase 
in  the  computational  load  if  all  the  angle  subintervals  are  processed  throughout.  However, 
it  has  been  found  that  generally  the  weighting  of  some  of  the  subintervals  rapidly  reduces 
to  zero.  In  such  cases,  the  corresponding  filters  can  be  removed  from  the  tracking  process 
without  loss  of  accuracy,  thereby  reducing  the  processing  requirement.  Thus,  a  weighting 
threshold  can  be  set  and  any  filter  corresponding  to  a  subinterval  with  a  weight  less  than 
the  threshold  may  be  removed  from  the  tracking  process. 

4.3  Regularised  Particle  filter 

The  particle  filter  (PF)  is  a  technique  for  implementing  a  recursive  Bayesian  filter  by 
Monte-Carlo  simulations  [6].  The  optimum  recursive  Bayesian  filter  in  the  Minimum  Mean 
Square  Error  (MMSE)  sense  is  the  mean  of  the  posterior  density  [13].  The  key  idea  of  the 
PF  is  to  represent  the  required  posterior  density  function  by  a  set  of  random  samples  or 
“particles”,  and  to  compute  the  state  estimate  based  on  these  samples.  As  the  number 
of  samples  becomes  very  large,  this  Monte-Carlo  characterisation  becomes  an  equivalent 
representation  to  the  usual  functional  description  of  the  posterior  pdf,  and  the  particle 
filter  estimate  approaches  the  optimal  Bayesian  estimate. 

4.3.1  SIR  algorithm 

In  order  to  recursively  compute  the  state  estimates  with  the  PF,  the  Monte  Carlo 
representation  of  p(sfc|Zfc)  has  to  be  propagated  in  time.  The  original  particle  (a.k.a. 
bootstrap)  filter,  proposed  by  Gordon  et  al.  [11],  was  using  the  Sampling  Importance 
Resampling  (SIR)  algorithm  to  propagate  and  update  the  particles.  Assume  at  time  tk— i 
we  have  a  set  of  independent  and  identically  distributed  (i.i.d)  random  samples 
from  the  posterior  pdf  p(sfc_i|Zfc_i).  Then  for  large  N  this  pdf  can  be  approximated  as: 

N 

p{sk-i  |Zfc_i)  «  6(sk-i  -  s|_!),  (77) 

i= 1 

where  <$(■)  is  the  Dirac  delta  function.  The  SIR  algorithm  propagates  the  random  sample 
{si-i)iLi  drawn  from  P(s*-ilz*-i),  to  a  new  random  sample  {sj.}^  which  is  approxi¬ 
mately  distributed  according  to  p(sfc|Zj,).  This  is  performed  as  follows.  According  to  the 
Bayes  rule  the  posterior  pdf  at  k  can  be  written  as  [13,  Sec. 6. 6]: 

p(s*,|Zfc)  <x  p(zk\sk) piskl'Zk-i)  (78) 

where  the  prediction  density  p(sfc|Zfc_i)  is  accomplished  via 

p(sfc|Z*_1)  =  J  p{sk\sk-i)p{sk-i\Zk-i)dsk-i.  (79) 

Using  (77)  in  (79),  the  prediction  density  can  be  approximated  as 

p(sjfc|Zfc_1)  «  ±  2>(sft|4-i).  (80) 
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The  sample  from  p(sjc\Zfc)  is  generated  in  two  stages,  namely  prediction  and  update. 
In  the  prediction  stage,  N  random  samples  are  generated  from  the  prediction  density 
p(sfc|Zfc_i)  given  in  (80).  This  is  done  by  passing  the  sample  {sj^}^  through  the  state 
dynamic  model  of  eq.(14).  Note  that  in  the  absence  of  process  noise  in  the  state  dynamics 
(the  case  we  are  considering  in  this  report)  the  propagation  of  particles  is  purely  deter¬ 
ministic.  The  resulting  random  sample  {s ls  approximately  distributed  according  to 
p(s}c\Zk~1).  In  the  update  stage,  the  new  sample  {sj^}^  is  generated  by  resampling  N 
times  from  an  approximate  discrete  representation  of  p( Sjfc|Z*.)  given  by: 


N 

P(sk  izfc)«£gjU(Sfc-sr) 

2=1 


(81) 


where 

a  =  pfaK) 
k  £f=1PW)' 


(82) 


In  other  words,  samples  are  generated  from  the  discrete  set  {s^}^  such  that 

Pr( sj.  =  s^)  =  q\.  Note  that  in  the  update  stage,  the  measurement  z*  is  used  to  evaluate 
the  likelihood  functions  p(z*.|s£l)  in  (82). 

A  common  problem,  particularly  prevalent  in  dynamical  systems  with  small  or  zero 
process  noise,  is  the  sample  impoverishment  phenomenon,  where  the  number  of  distinct 
particles  at  each  time  step  monotonically  decreases.  If  the  problem  is  not  addressed 
properly,  it  will  lead  to  “particle  collapse” ,  where  all  N  particles  occupy  the  same  point 
in  the  state  space,  giving  a  poor  representation  of  the  posterior  density. 


4.3.2  Regular  isat  ion 

An  effective  scheme  that  overcomes  the  sample  impoverishment  problem  is  known  as 
the  regularisation  of  the  particle  filter  [16].  The  resulting  filter,  the  Regularised  Particle 
filter  (RPF),  differs  from  the  SIR  in  the  resampling  step  only.  The  RPF  resamples  from  a 
continuous  approximation  of  the  posterior  density  p(sjt)Z^),  while  the  SIR  resamples  from 
the  discrete  approximation  (81).  Specifically,  in  the  RPF,  samples  are  drawn  from  the 
approximation: 

N 

p(sfc|Zfc)«^4^(sfc-Sn  (83) 

2  =  1 

where 

**<s> = ikK$  m 

is  the  re-scaled  kernel  density  RT(*),  h  >  0  is  the  kernel  bandwidth  (a  scalar  parameter), 
and  ns  is  the  dimension  of  the  state  vector  s.  The  kernel  density  is  a  symmetric  probability 
density  function  such  that 
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The  kernel  K(-)  and  bandwidth  h  are  chosen  so  as  to  minimise  the  mean  integrated 
square  error  between  the  true  posterior  density  and  the  corresponding  regularised  empirical 
representation  in  (83).  In  the  special  case  of  an  equally  weighted  sample,  i.e.  with  q\  = 
l/N  for  i  =  1,2,. ..,7V,  the  optimal  choice  of  the  kernel  is  the  Epanechnikov  kernel  [16], 


10 


if||s||<l 

otherwise. 


(86) 


When  the  underlying  density  is  Gaussian  with  a  unit  covariance  matrix,  the  optimal  choice 
for  the  bandwidth  is 

hopt  =  A(k)N~^r*  (87) 

A{k)  =  [8c-1(na  +  4)(2^)ns];^3  (88) 

where  Cn,  is  the  volume  of  the  unit  hypersphere  in  R”» .  The  general  formula  for  Cn,  is 
given  in  [3,  p.96].  For  the  state  vector  as  in  (3),  na  =  4  and  C4  =  n2/2. 

The  implementation  of  the  RPF  is  described  next.  The  prediction  step  is  the  same  as 
in  the  SIR  algorithm,  and  it  creates  a  random  sample  {s^}^.  Now  instead  of  resampling 
directly  from  {s ^}^=1  as  done  by  SIR,  the  following  algorithm  is  implemented: 


1.  Draw  N  indices  from  a  set  so  that  Pr{h  =  j}  =  q{,  where 

2.  Calculate  the  empirical  covariance  matrix  S(.  of 

3.  For  each  i  =  1, . . .  ,N 

a.  Generate  a  random  variable  e*  ~  K  from  Epanechnikov  kernel. 

b.  Compute 

si  =  sf  +  hopeDfce'  (89) 

where  is  the  square  root  of  the  empirical  covariance  matrix  S£,  i.e.  DfcDj  = 

S£- 

By  following  this  procedure  we  generate  a  random  sample  {s[.}iLi  drawn  from  (83). 


5  Algorithm  Performance  and  Comparison 

The  statistical  performance  of  the  TMA  algorithms  is  analysed  by  Monte  Carlo  simu¬ 
lations.  The  trajectories,  the  measurement  covariances  R  and  Ro,  the  sampling  interval 
Ts  and  the  initial  delay  Td  have  the  same  values  as  described  in  Section  3.2.  The  resulting 
error  curves  are  compared  to  the  theoretically  derived  Cramer-Rao  lower  bounds. 

Define  the  estimation  error  as: 

ek  =  ^k\k  ~~  Sk-  (90) 
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The  performance  of  the  TMA  algorithms  is  measured  by  the  root  mean  square  (RMS) 
error.  For  component  j  of  the  state  vector,  the  RMS  error  is  defined  as 


a[ek]j  ~~  \/®  {[efc]j  E{[e^]j}} 

and  approximated  with  M  Monte  Carlo  runs  as  follows: 


a[e*]i 


N 


i  M  (  i  M 

fc=l  l.  k- 1 


(91) 


(92) 


The  theoretical  CR  bound  is  derived  earlier  for  unbiased  estimators.  Its  square  root 
represents  the  lower  bound  for  RMS  error  of  any  unbiased  TMA  algorithm,  i.e. 


a[et].  >  y'CRLBds^}  (93) 

where  CRLB^sj^*}  is  given  by  eq.(53).  In  all  simulation  results  we  used  M  =  100  inde¬ 
pendent  Monte  Carlo  runs. 


5.1  Performance  of  the  MLE 

The  statistical  performance  of  the  Maximum  likelihood  estimator  is  shown  in  Figure  6 
(x  and  y  position  only,  i.e.  j  =  1, 3  in  (92)).  Since  the  MLE  is  not  a  Bayesian  estimator,  its 
performance  is  compared  to  the  CRLB  with  no  use  of  prior  knowledge.  Figure  6  confirms 
the  theoretical  assertion  that  the  MLE  for  range-only  TMA  is  an  efficient  estimator:  it 
has  no  bias  and  its  RMS  error  attains  the  square-root  of  the  calculated  CRLB. 


Figure  6:  The  statistical  performance  of  the  MLE  against  the  CRLB:  (a)  x  position  and  (b)  y 
position 


The  main  disadvantage  of  the  MLE  is  its  computational  complexity.  Not  only  is  the 
numerical  minimisation  of  the  log- likelihood  function  in  (61)  a  computationally  demanding 
operation,  but  also  at  every  time  index  k  the  MLE  must  operate  on  the  accumulated  set  of 
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all  previous  measurements4  Zk.  Using  the  same  scenario  as  before,  we  compared  the  CPU 
time  required  by  the  MLE  algorithm  to  that  of  the  EKF  algorithm.  Both  algorithms  were 
implemented  in  MATLAB®,  and  the  comparison  is  performed  at  time  index  k  =  240.  It 
appears  that  the  MLE  requires  580  times  more  CPU  time  than  the  EKF. 


5.2  Performance  of  the  EKF 

The  statistical  performance  of  the  EKF  based  range-only  TMA  algorithm  is  shown  in 
Figure  7  (x  and  y  position  only,  i.e.  j  -  1,3  in  (92)).  The  EKF  is  a  Bayesian  estimator, 
hence  its  performance  is  compared  to  the  CRLB  with  the  prior  knowledge  of  ov  =  10  m/s. 
For  comparison  sake,  the  MLE  performance  curves  are  indicated  by  a  dotted  line. 


Figure  7:  The  statistical  performance  of  the  EKF  against  the  CRLB  (with  prior  knowledge  of 
uv  =  10  m/s):  (a)  x  position  and  (b)  y  position.  The  dotted  line  is  the  MLE  error  curve. 


Observe  from  Figure  7  that  the  EKF  error  is  initially  smaller  than  the  MLE  error, 
since  the  EKF  uses  the  prior  knowledge  about  the  maximum  target  velocity.  As  the  time 
progresses,  however,  the  EKF  quickly  departs  from  the  CRLB,  and  in  general  demonstrates 
a  slower  convergence  than  the  MLE.  This  is  particularly  pronounced  in  Figure  7.(b).  The 
suboptimal  performance  of  the  EKF  is  due  to  the  approximation  (linearisation)  of  the 
non-linear  measurement  equation. 

Angle  parametrisation  described  in  Section  4.2.2  effectively  resolves  this  problem.  The 
results,  obtained  using  Nf  =  5  EKF  filters,  are  shown  in  Figure  8.  Again,  the  performance 
of  the  AP-EKF  is  compared  to  the  CRLB  with  the  prior  knowledge  of  av  =  10  m/s,  and 
the  MLE  error  curves  are  indicated  by  a  dotted  line  for  comparison.  Initially,  the  AP-EKF 
outperforms  the  MLE  (by  using  the  prior  knowledge  of  the  maximum  target  velocity),  and 
subsequently  its  performance  is  very  close  to  that  of  the  MLE.  Overall,  this  represents  a 
huge  improvement  over  the  EKF.  Our  implementation  of  the  AP-EKF  (in  the  worst  case) 
required  Nf  =  5  times  more  CPU  time  than  the  EKF  and  580/A/  =  116  times  less  CPU 
time  than  the  MLE. 

4 By  restricting  the  MLE  to  operate  on  the  sliding  window  of  previous  measurements,  the  efficiency  of 
the  MLE  is  lost. 
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Figure  8:  The  statistical  performance  of  the  AP-EKF  against  the  CRLB  (with  prior  knowledge 
of  av  —  10  m/s):  (a)  x  position  and  (b)  y  position.  The  dotted  line  is  the  MLE  error  curve. 

It  should  be  emphasised  that  in  all  simulations  the  value  of  ere  =  \/[Rol22  =  3° 
was  used.  However,  if  the  accuracy  of  the  initial  angular  measurement  is  higher  (i.e. 
a#  <  1°),  the  EKF  based  TMA  actually  attains  its  corresponding  CRLB  and  the  angle 
parametrisation  of  the  EKF  is  unnecessary. 

5.3  Performance  of  the  RPF 

The  Regularised  Particle  filter  was  included  for  consideration  in  this  study  of  range- 
only  TMA  algorithms,  mainly  because  it  approaches  in  theory  (with  large  number  of 
particles  N)  the  optimal  recursive  Bayesian  estimator  for  what  is  essentially  a  non-linear 
dynamic  estimation  problem. 

In  order  to  be  fair,  the  RPF  should  be  compared  to  the  EKF  filter  for  range-only 
TMA,  because  the  angle  parametrisation  technique  discussed  above  in  the  context  of  EKF 
is  equally  applicable  to  the  RPF.  Figure  9  shows  the  statistical  error  performance  of  the 
RPF  based  TMA  algorithm  ( x  and  y  positions  only)  with  N  =  2000,  5000  and  10000 
particles.  For  comparison  sake,  the  CRLBs  (using  the  prior  knowledge  of  av  =  10  m/s) 
and  the  EKF  error  curves  are  shown  as  well.  Observe  that  as  the  number  of  particles 
N  is  increased,  the  performance  of  the  RPF  improves,  and  in  the  limit  approaches  the 
CRLBs.  This  improvement,  however,  is  at  the  expense  of  a  corresponding  increase  in 
the  computational  load  of  the  RPF.  Thus  the  RPF  with  2000,  5000  and  10000  particles 
requires  30,  70  and  134  times  more  CPU  time  respectively,  than  the  EKF.  Note  also  that 
at  the  end  of  the  observation  period,  the  RMS  error  of  the  EKF  is  slightly  better  than  the 
RMS  error  of  the  RPF.  This  is  due  to  the  sample  impoverishment  in  the  particle  filter, 
which  has  an  effect  towards  the  end  of  the  observation  period  despite  regularisation. 

Finally,  we  make  an  overall  comment  with  respect  to  the  performance  of  various  range- 
only  TMA  algorithms.  It  appears  that  the  AP-EKF  would  be  the  strongest  candidate  for 
the  considered  application.  The  AP-EKF  combines  the  satisfactory  statistical  perfor¬ 
mance,  robustness  to  the  initial  measurement  accuracy  and  the  speed  of  execution.  An- 
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Figure  9:  The  statistical  performance  of  the  RPF  with  N  =  2000,  5000  and  10000  particles. 
The  CRLBs  (using  prior  knowledge  of  av  =  10  m/s)  are  shown  with  dashed  lines;  the  EKF 
curves  are  plotted  with  dotted  lines,  (a)  x  position  and  (b)  y  position. 


other  advantage  of  the  AP-EKF  over  the  MLE  is  that  it  can  directly  incorporate  a  small 
level  of  process  noise  in  the  dynamic  equation  (14).  This  can  be  important  for  modelling 
small  changes  in  the  target  velocity  resulting  from  the  rough  sea  in  windy  conditions. 


6  Application  to  Ingara  ISAR  Data 

6.1  The  ISAR  mode  of  Ingara  radar  for  data  collection 

Inverse  Synthetic  Aperture  Radar  (ISAR)  is  a  technique  for  generating  high-resolution 
radar  images  in  range  and  Doppler  (cross-range).  This  is  achieved  by  sampling  radar  echoes 
from  a  target  and  coherently  processing  blocks  of  echoes  to  form  ISAR  images.  The  high 
resolution  in  range  is  achieved  by  using  a  stretch  waveform  while  the  high  resolution  in 
cross-range  is  achieved  by  processing  a  large  aperture  of  data.  A  comprehensive  treatment 
of  ISAR  theory  is  given  in  [21]. 

The  ISAR  mode  on  Ingara  has  been  developed  by  Surveillance  Systems  Division  of 
DSTO.  A  series  of  trials  have  been  conducted  to  build  up  a  database  of  ISAR  imagery 
of  a  representative  collection  of  target  types  to  assist  in  the  process  of  operator  training 
and  for  further  research  into  automatic  target  recognition.  The  ISAR  mode  of  the  Ingara 
radar  collects  and  coherently  processes  high  resolution  range  profiles  into  ISAR  images  for 
real-time  display.  This  is  achieved  by  “spotlighting”  a  maritime  surface  vessel  during  the 
data  collection.  The  typical  radar  parametres  used  for  the  ISAR  mode  axe  shown  in  Table 
1. 

The  application  of  range-only  tracking  algorithms  to  a  set  of  real  ISAR  data  will  be 
described  below.  The  Ingara  radar  is  installed  on  a  Beechcraft  KingAir  350  aircraft,  and 
targets  of  opportunity  were  used  in  data  collection. 
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Table  1:  Parameters  used  in  Ingara  ISAR  mode 


Radar  Parameters 

A/D  sampling  rate 

100  MHz 

RF  Bandwidth 

400  MHz 

Pulse  Width 

18  fis 

Number  of  samples 

2048 

Slant  range  resolution 

0.33  m 

Typical  target  range 

20  km 

PRF 

600  Hz 

Transmitted  power 

5  kW 

Polarisation 

VV 

Rut  first,  we  review  the  coordinate  conversion  techniques,  necessary  for  transformation 
of  radar  measurements  to  the  local  tangential  plane  Cartesian  coordinate  system,  which 
was  assumed  from  Section  2  to  Section  5. 


6.2  Coordinate  transformations 

Four  coordinate  systems  are  used  in  the  analysis  of  ISAR  real  data  (Fig.  10).  The  own- 
ship  coordinates,  supplied  by  the  on-board  global  positioning  system  (GPS),  are  given 
in  the  geodetic  coordinate  system.  The  TMA  problem  and  the  corresponding  algorithms 
(Sections  2-5)  are  defined  in  the  tangential- plane  Cartesian  coordinate  system.  The  radar 
measurements  (range,  range-rate,  initial  azimuth)  are  collected  in  the  local  polar  coordi¬ 
nates.  Finally  the  transformation  between  the  geodetic  and  the  tangential- plane  Cartesian 
coordinate  systems  is  done  via  Earth- Centered  Earth-Fixed  (ECEF)  coordinates.  The  con¬ 
versions  between  the  coordinates  is  described  below  following  [9,  Chap. 2].  Note  however 
that  the  conversion  from  the  local  polar  to  the  local  Cartesian  coordinates  is  not  necessary 
as  the  TMA  algorithms  use  the  radar  measurements  directly. 


Figure  10:  Coordinate  systems. 


Geodetic  to  ECEF  and  vice  versa.  In  the  geodetic  coordinate  system,  a  point  on 
the  ellipsoidal  earth  is  represented  by  a  triple  (A,  </>,  /i),  where  A  is  latitude,  <f>  is  longitude 
and  h  is  the  altitude  above  the  reference  ellipsoid.  In  the  ECEF  rectangular  coordinates 
(X,  y,  Z),  the  origin  is  in  the  center  of  the  earth.  The  X  axis  extends  through  the  inter- 
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section  of  the  prime  (0°  longitude)  meridian  and  the  equator  (0°  latitude).  The  Z  axis 
extends  through  the  north  pole  (coincides  with  the  earth’s  spin  axis).  The  Y  axis  com¬ 
pletes  the  right-handed  coordinate  system,  passing  through  the  equator  and  90°  longitude. 
The  earth  geoid  is  approximated  by  an  ellipsoid  of  revolution  about  its  minor  axis.  The 
parameters  of  the  chosen  approximating  ellipsoid  must  be  defined  for  a  geodetic  system, 
and  in  our  case  this  is  the  WGS-84  ellipsoid.  Its  defining  parameters  are  the  semimajor 
and  semiminor  axis  lengths  a  and  b  respectively  [9,  p.26].  The  coordinate  transformations 
from  geodetic  to  ECEF  coordinates  are  as  follows: 

X  =  (N  +  h)  cos  A  cos  4> 

Y  —  (N  +  h)  cos  A  sin  <f>  (94) 

Z  -  [AT(1  -  e2)  +  h]  sin  A  (95) 


where  e  is  the  eccentricity  of  the  ellipsoid  (for  WGS-84,  e  =  0.0818)  and 


N(  A)  = 


a 

\/l  —  e2  sin2  A 


(96) 


The  conversion  from  the  ECEF  to  the  geodetic  coordinates  is  more  complicated  and  we 
used  the  iterative  solution  described  in  [9,  p.28]. 


ECEF  to  tangent-plane  and  vice  versa.  The  local  tangent-plane  coordinate  system 
is  the  east-north-up  rectangular  coordinate  system  we  often  refer  to  in  our  everyday  life. 
It  is  determined  by  the  fitting  of  a  tangent  plane  to  the  surface  of  the  earth  at  some  conve¬ 
nient  (reference)  point  for  local  measurements  (typically  in  vicinity  of  the  data  collection 
scenario).  This  reference  point  is  the  origin  of  the  local  frame.  The  x  axis  points  to  true 
east;  the  y  axis  points  north  and  the  2  axis  completes  the  right-handed  coordinate  system 
pointing  up.  For  a  moving  own-ship,  the  tangent-plane  origin  is  fixed. 

The  transformation  matrix  for  transformation  from  ECEF  to  tangent-plane  represen¬ 
tation  is  given  by  [9,  p.35]: 


C  = 


—  sin  cp 
—  sin  A  cos(j> 
cos  A  cos  0 


cos  <f>  0 

—  sin  A  sin  cj>  cos  A 
cos  A  sin  4>  sin  A 


(97) 


If  the  coordinates  of  the  reference  point  in  ECEF  are  (X0,  To,  Z0),  then  the  transformation 
from  an  ECEF  point  (X,  Y,  Z)  to  the  local  tangent-plane  coordinates  ( x ,  y,  z)  is  as  follows: 


x 

y 

z 


The  inverse  transformation  is  then: 


'  X  ■ 
Y 

_ 

'  *o  ' 
To 

+  CT 

X 

y 

Z 

Zo 

z 

(98) 


(99) 
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6.3  Pre-processing  of  the  ISAR  data 

The  target  measurements  which  feed  the  TMA  algorithms  are  obtained  by  pre-processing 
a  large  number  of  range  profiles  as  follows.  The  data  is  divided  into  contiguous  blocks 
corresponding  to  the  desired  update  rate  for  the  tracker.  The  estimates  of  target  range 
and  range-rate  (referred  to  earlier  as  measurements  r &  and  r&)  are  computed  for  each 
block.  In  this  particular  case  we  adopt  the  sampling  rate  of  1  Hz  (i.d.  Ts  —  1  s)  and  since 
the  radar  PRF  is  600  Hz,  a  block  size  of  600  range  profiles  is  used  to  compute  each  pair 
rk  and  ?V  An  initial  range  estimate  is  formed  by  averaging  the  first  few  range  profiles 
at  the  start  of  the  very  first  block  and  by  computing  its  “centroid”  (a  weighted  mean). 
This  is  only  done  for  the  first  block  and  represents  the  nominal  target  range  at  the  start 
of  the  ISAR  data.  Each  subsequent  range  profile  (from  the  second  onwards)  is  correlated 
with  the  first  profile  and  the  correlation  lag  (the  shift  in  the  position  of  the  maximum 
in  the  correlation  sequence)  is  due  to  target  motion.  Only  shifts  of  ±1  range  bin  from 
one  profile  to  the  next  are  allowed  because  the  target  cannot  physically  move  faster.  Any 
shift  of  more  than  ±1  bin  is  replaced  by  a  zero  shift  because  this  indicates  bad  or  missing 
data  in  that  profile  (e.g.  due  to  interference  or  hardware  fault).  Having  computed  the 
target  shift  for  each  profile  by  correlation,  we  fit  through  them  a  line  by  the  least-squares 
method.  The  range-rate  r*  is  calculated  by  taking  the  slope  of  that  line.  The  target  range 
from  the  middle  profile  of  the  block  is  used  as  the  range  estimate  r^.  For  each  subsequent 
block  the  same  procedure  is  repeated,  except  that  the  initial  range  estimate  of  the  block 
is  not  formed  by  the  centroid  of  the  averaged  range  profiles,  but  simply  by  using  the  range 
estimate  from  the  last  range  profile  of  the  previous  block. 

With  the  described  procedure  for  calculation  of  r &  and  r&,  the  assumptions  made  in 
relation  to  the  measurement  noise  wjt  in  (6),  become  satisfied  in  a  broad  sense.  In  addition, 
the  measurements  of  range  r ^  and  range-rate  r k  are  uncorrelated  (i.e.  the  covariance  matrix 
R  is  diagonal). 


6.4  Range-only  TMA  from  the  ISAR  data 

The  geographic  location  of  the  data  collection  scenario  is  shown  in  Figure  11.  The 
initial  observer  position  and  the  initial  target  location  measurement  are  indicated  by  0 
and  ★  respectively.  The  reference  point  (the  origin  of  the  local  coordinates)  is  shown  as  □. 
The  observer  trajectory  after  the  time  delay  of  T =  65  sec  is  plotted  by  a  solid  line.  The 
sampling  interval  (i.e.  the  time  increment  of  the  sliding  window  used  for  pre-processing 
range  profiles)  is  Ts  =  1  sec.  The  observer  height  was  constant  at  210  m. 

The  own-ship  position  and  velocity,  computed  from  the  supplied  INS/GPS  data  and 
transformed  into  the  local  tangent  plane,  are  shown  in  Figure  12.  Finally  the  range 
and  the  range  rate  measurements  are  displayed  in  Figure  13.  Occasional  outliers  in  the 
measurements  are  filtered  out  using  a  median  filter  of  windowlength  3.  Since  the  target 
range  is  approximately  20  km  (see  Figure  13. a),  in  all  our  TMA  calculations  we  neglect 
the  own-ship  height,  and  assume  that  the  observer  and  the  target  are  co-planar.  The  error 
in  the  measured  range  due  to  this  approximation  is  less  than  1  m. 

The  ISAR  data  were  collected  using  the  targets  of  opportunity.  While  the  true  target 
state  over  the  observation  period  is  unknown,  we  have  two  means  of  testing  and  verifying 
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<0  Ownship  Init  Posit 
——  Ownship  Trajectory 
☆  Target  Init  Posit 


Figure  11:  The  geographic  location  of  the  data  collection  scenario 


the  TMA  algorithm  performance. 

First  we  apply  the  MLE  TMA  algorithm  on  the  entire  measurement  set  Z k,  (where  K 
is  the  total  number  of  available  measurements)  and  then  compare  the  measured  range  and 
range-rate  with  what  would  have  been  measured,  had  the  estimated  target  state  been  cor¬ 
rect.  The  result  of  this  test  with  R  =  diag[(20m)2;  (3m/s)2]  and  Ro  =  diag[(20m)2;  (1°)2] 
is  shown  in  Figure  14.  We  note  in  the  range  plot  (Figure  14.a)  that  the  two  curves  appear 
to  be  identical,  except  at  to  =  0,  where  the  discrepancy  is  about  200  m.  This  is  due  to  the 
lack  of  own-ship  motion  compensation  during  the  period  of  time  Td  required  for  the  radar 
mode  switching  (the  own-ship  data  were  unavailable  during  this  interval).  The  range-rate 
curves  in  Figure  14.  )  are  also  very  close,  well  within  the  assumed  standard  deviation  of 
the  range-rate  error  of  \/[R]22  =  3  m/s.  This  result  gives  us  confidence  that  the  MLE 
estimate  taken  using  the  entire  measurement  set  Z k  is  accurate. 

The  second  test  uses  the  information  about  the  true  target  heading,  which  was  mea¬ 
sured  independently  as  287°  clockwise  positive  with  respect  to  North.  Assuming  the  target 
position  measurement  at  to  (that  is  zo)  is  exact  and  using  the  estimate  of  the  target  speed 
(5.9  m/s)  obtained  from  the  MLE  state  estimate,  we  obtain  what  we  refer  to  as  the  in¬ 
complete  truth  data  about  the  target.  The  target  trajectory  based  on  incomplete  truth 
is  shown  in  Figure  15  by  a  dotted  line.  For  comparison,  Figure  15.a  further  displays  the 
estimated  target  trajectories  of  the  progressive  MLE  (applied  to  Zfe,  k  =  1, . . . ,  K);  the 
angle-parametrised  EKF  and  the  Regularised  PF.  The  AP-EKF  was  using  Nf  =  5  filters 
and  both  AP-EKF  and  RPF  used  <jv  =  10  m/s.  Note  that  all  three  algorithms  fairly  accu¬ 
rately  follow  the  indicated  target  heading.  The  MLE  and  the  AP-EKF  converge  towards 
the  true  trajectory  after  approximately  2.5  minutes,  while  the  RPF  convergence  is  some¬ 
what  slower.  For  completeness,  Figure  15.b  displays  the  MLE  estimated  target  trajectory 
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Figure  12:  Own-ship  position  and  velocity  in  the  local  coordinates 


on  a  geographic  map  (in  geodetic  coordinates).  Finally  Figure  15. c  repeats  the  estimated 
target  trajectory  obtained  with  the  AP-EKF,  but  it  incorporates  the  2-sigma  uncertainty 
ellipses.  Note  that  these  ellipses:  (i)  always  include  the  true  target  position;  (ii)  reduce 
in  size  as  time  progresses;  (iii)  indicate  a  large  uncertainty  in  angle  (not  measured)  but 
small  in  range;  (iv)  have  a  long  axes  always  orthogonal  to  the  “line-of-sight” . 


7  Summary 

The  report  describes  an  invetsigation  of  the  problem  of  target  motion  analysis  from 
range  and  range-rate  measurements.  The  problem  is  of  importance  for  tracking  surface 
vessels  while  collecting  data  in  the  ISAR  mode  of  the  Ingara  multi-mode  radar.  The  re¬ 
port  derived  the  bounds  of  performance  and  considered  three  algorithms  for  target  motion 
analysis.  The  proposed  algorithms  include  the  Maximum  Likelihood  estimator,  the  Ex¬ 
tended  Kalman  filter  with  and  without  angle  parametrisation  and  the  Regularised  Particle 
filter.  Considering  both  the  statistical  and  computational  performance  of  the  proposed 
algorithms,  this  study  recommends  the  angle-parametrised  EKF  as  the  preferred  choice 
for  implementation  in  an  operational  system. 

The  discussed  TMA  algorithms  have  been  applied  successfully  to  one  set  of  real  ISAR 
data  collected  in  the  recent  trials  with  the  Ingara  radar.  The  results  indicate  that  target 
motion  analysis  based  on  range  and  range-rate  measurements  can  converge  towards  a 
steady  state  in  less  than  three  minutes  for  a  typical  scenario  of  interest  for  ISAR  data 
collection. 

Future  work  would  involve  real-time  implementation  of  the  recommended  range-only 
tracking  algorithm  on  Ingara  system. 
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Figure  13:  Target  range  and  range-rate  measurements 
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contributions.  First,  the  theoretical  Cramer-Rao  bound  for  the  performance  of  an  unbiased  range-only 
tracking  algorithm  is  derived.  Second,  three  algorithms  for  target  motion  analysis  (using  range  and 
range-rate  measurements  only)  are  developed  and  compared  to  the  theoretical  bounds  of  performance. 
The  three  algorithms  are:  the  Maximum  Likelihood  estimator,  the  Extended  Kalman  filter  and  the 
Regularised  Particle  filter.  Finally,  the  report  presents  the  results  of  the  application  of  the  developed 

theory  to  the  ISAR  data  collected  in  the  recent  trials  with  the  Ingara  radar. _ 
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